%dXt = sigma(t)^2dt 
%sigma(t)=OU
clear all
% index J for X
betaK=15; IIpart=sparse(index2(betaK));
% time T
T = 1/2;
% time mesh 
num=50; tspan = linspace(0,T,num+1);
%sigma
s=sigmahandle; a=1;b=1;m=1;sig0=1;
%options = odeset('AbsTol',1e-12,'RelTol',1e-6);
options = [];
row=1+2*betaK+betaK*(betaK-1)/2; X0=zeros(row,1); 
% solver of X
disp('solving X1')
tic
[~,X1] = ode45(@randucousqdt,tspan,X0,options,T,a,b,m,sig0,betaK,s,IIpart);

%dXt = sigma(t)dW 
%sigma(t)=OU

% index J for X
alphaK=20; %betaK=15; 
%sigma
sig=sigmahandle; %a=1;b=1;m=1;sig0=1;
s  =mihandle;
row=(betaK+1)*(alphaK+1); X0=zeros(row,1); 
% solver of X
disp('solving X2')
[time,X2] = ode45(@randucoudw,tspan,X0,options,T,a,b,m,sig0,alphaK,betaK,sig,s);


%%%the WCE price
price=zeros(51,1);
Npath=20000;
eta=randn(Npath,betaK);
etasquare=1/sqrt(2)*(eta.^2-1);
etaeta=etaietaj(eta);
xi=randn(Npath, alphaK);

for iter=1:51
% asset=zeros(Npath,1);
% for ii=2:betaK+1
%     asset=asset+X1(iter,ii)*eta(:,ii-1);
% end
% for ii=betaK+2:2*betaK+1
%     asset=asset+X1(iter,ii)*etasquare(:,ii-betaK-1);
% end
% for ii=2*betaK+2:2*betaK+1+betaK*(betaK-1)/2
%     asset=asset+X1(iter,ii)*etaeta(:,ii-2*betaK-1);
% end
% asset=asset+X1(iter,1);

asset=pricex1(eta,etasquare,etaeta,X1(iter,:));
% 
% asset=X2(iter,1)*ones(Npath,1);
% for jj=1:betaK
%     asset=asset+X2(iter,jj+1)*eta(:,jj);
% end
% for ii=1:alphaK
%     asset=asset+X2(iter,ii*(betaK+1)+1)*xi(:,ii);
%     for jj=1:betaK
%         asset=asset+X2(iter,ii*(betaK+1)+1+jj)*xi(:,ii).*eta(:,jj);
%     end
% end
asset=asset+pricex2(eta,xi,X2(iter,:));
price(iter)=sum(abs(asset))/Npath;
end
toc

clear eta etasquare etaeta asset xi X1 X2 X0 IIpart s sig tspan

disp('MC price')
tic
Npath=20000;

randseed =randn(5000,Npath);
sampleou=oupath(randseed,1,1,1,1);
randseed =randn(5000,Npath); % because of uncorrelated
sampleX =mcprice(sampleou,randseed);
tmp =sum(abs(sampleX),2)/Npath;
mtclprice=tmp(1:100:5001);
toc

% tic
% Npath=10000;
% sampleX=zeros(5001,Npath);
% for ii=1:5
% 
% 
% randseed =randn(5000,Npath);
% sampleou=oupath(randseed,1,1,1,1);
% randseed =randn(5000,Npath); % because of uncorrelated
% sampleX =sampleX+mcprice(sampleou,randseed);
% 
% end
% sampleX=sampleX/5/Npath;
% tmp =sum(abs(sampleX),2);
% mtclprice=tmp(1:100:5001);
% 
% 
% % tic
% % intsigmasq=zeros(51,1);
% % for pn=1:Npath
% %     samplepath=getapath(5001,'ou');
% %     sigmasq=samplepath.^2; 
% %     for jj=1:51
% %         intsigmasq(jj)=intsigmasq(jj)+midptrule(sigmasq(1:(jj-1)*100+1),1/10000);
% %     end
% % end
% % mcpart=zeros(51,1);
% % for pn=1:Npath
% %     samplepath=getapath(5001,'sigmadW');
% %     sigmadW=abs(samplepath); 
% %     for jj=1:51
% %         mcpart(jj)=mcpart(jj)+sigmadW((jj-1)*100+1);
% %     end
% % end
% % mcpart=mcpart/Npath;
% % 
% % intsigmasq=intsigmasq/Npath;
% toc

hold on
plot(time, price);
plot(time, mtclprice,'r')



